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Abstract 


An algorithm is presented which solves the multidimensional diffusion equation 
on complex shape* to 4 tk -order accuracy sad is asymptotically stable in time. This 
bounded-error result is achieved by constructing, cm a rectangular grid, a differentiation 
matrix whose symmetric part is negative definite. The differentiation matrix accounts 
for the Dirichlet boundary condition by imposing penalty Uhe terms. 

Numerical examples in 2-D show that the method is effective oven where standard 
schemes, stable by traditional definitions, foO. 
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1 Introduction 


Recently there has been renewed interest in finite- difference algorithms of high order of 
accuracy (4* and above), both for hyperbolic and parabolic p.d.e’s (see for example, [1], [2], 
[3] ). The advantages of high-order accuracy schemes, especially for truly time dependent 
problems, are often offset by the difficulty of imposing stable boundary conditions. Even 
when the scheme is shown to be G.K.$.-stable the error may increase exponentially in time. 

This paper is concerned with 4 u *-order approximations to the long time solutions of the 
diffusion equation in one and two dimensions, on irregular domains. By an irregular domain, 
we mean a body whose boundary points do not xnntide with nodes of a rectangular mesh. 

In section 2 we develop the theory for the one-dimensional semi-discrete system resulting 
from the spatial differentiation used in the finite difference algorithm. Energy methods are 
used in conjunction with “SAT 1 type terms (see [1]), in order to find boundary conditions 
that preserve the accuracy of the scheme while constraining an energy norm of the error to 
be temporally bounded for all t > 0 by a constant proportional to the truncation error. 

In section 3 it is shown bow the methodology developed in section 2 is used as a building 
block for the multi-dimensional algorithm, even for irregular shapes containing "holes.” 
Section 4 presents numerical results in two space dimensions illustrating the long-time 
temporal stability of the method, in contradistinction to "standard” methods for cartesian 
grid on irregular shapes. 
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2 The One Dimensional Case 


We consider the Mowing problem 
dll 

di~ k lh 5 +/ (*’ , ) ? r tS*<ri», f>0, k > 0 

tt(* t 0) * «o(*) 

«(T L ,t) = 9L(i) 

* $*(<) 

and /(*,<) € C 4 . 

Let us spatially discritixe (2.1a) on the Mowing uniform grid: 



*1 *i *S X J-I *1 *J*t ^ 

Figure 1: One dimensional grid. 


Note that the boundary points do not necessarily coincide with x» and xs- Set * J+1 -Zj = A, 

— 1 ? o<*>t<i; Tr — */t — mb, o^7a<i. 

The projection unto the above grid of the exact solution u(x,t) to (2.1), is u,(t) » 
= **(<). Let D be a matrix representing the second partial derivative with respect to 
x, at internal points without specifying yet how it is being buih. Then we may write 

!«i(l)-MM<) + B + T) + f(<) (22) 


( 2 . 1 ») 

(2.1b) 

(2.1c) 

(2* Id) 
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where T is the truncation error due to the numerical differentiation and f(t) = f(xj,t), 
1 < J < iV. The boundary rector B has entries whose values depend on 9l,9r, 1l>1r in 
such a way that Du + B represents the derivative everywhere to the desired accuracy. 
The standard way of finding a numerical approximate solution to (2.1) is to omit T from 
(2.2) and solve 

|v(t)-*(0v(t) + B) + f(() (2.3) 

where v(t) is the numerical approximation to the projection u(t). An equation for the 
solution error vector, e{<) - u(t) - v(t), can be found by subtracting (2.3) from (2.2): 

jf-ifcW + HW (2.4) 

Our requirement for temporal stability is that || e |, tbe L, norm of e, be bounded by a 
“constant* proportional to k m (m being the spatial order of accuracy) for all t < oo. Note 
that this definition is more severe than either the G.K.S. stability criterion (4) or the definition 
in [1]. 

It can be shown that if D is constructed in a standard manner, i.e., the numerical second 
derivative is symmetric away from the boundaries, and near the boundaries one uses non 
symmetric differentiation, then there are ranges of values of tr and 74 for which D is 
not negative definite. Since in the multi*dimensional case one may MWMintff ill values of 
0 < 7t>7« < 1? this is unacceptable. 

The rest of this section is devoted to the construction of a scheme of 4* order spatial 
accuracy, which is temporally stable for all 71,7*. 
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The basic idea is to use at penalty-like tens as in the SAT procedure of ref [ 1 ]; here, 
however, it will be modified and applied in a different manner. 

Note first that the solution projection Uj(i) satisfies, besides (2.2), the following differ- 
ential equation: 

ij-U>u + iT, + f(l) (2.5) 

at 

where now D is indeed a differentiation matrix, that does not use the boundary values, and 
therefore T t ^ T but it too is a truncation error due to differentiation. 

Next let the semi-discrete problem for v(t) be, instead of (2.3), 

i * k[Dv - tl(Alv - gt) - t*(A*y - git)) + f(t) (2.6) 

where gt, * ( 1 , . . . , l) T pt(t); gjt * ( 1 , . . . * 1 ) T pj*(f ). are vectors created from the left and 
right boundary values as shown. Hie matrices At and Ar are defined by the relations: 

Ain m gt - T*; * g* - T*, (2.7) 

i.e., each row in /U,(A R ) is composed of the coefficients extrapolating u to its boundary value 
gt(gfl), at r t (r R ) to within the desired order of accuracy. (The error is then Ti,(T R ).) The 
diagonal matrices rj, and tr are given by 

ti m dSag(n* 9 «lt •»•«**»); r R * diag(r* t ,...r* K ) ( 2 . 8 ) 

Subtracting (2.6) from ( 2 . 5 ) we get 

5 = k\Df- 1 iA L i- T H A H i+ T,] 

4 


(2.9) 


I 


where 

Tt = T, + tlT l + r*T* 

I 

l Taking the scalar product of ? with (2.9) one gets: 

= *(«>(# - n.A L - r«.4«)e) + *(? t T,) 

I = K?,MQ + k& Tj) 

I 

j We notice that (?, MS) is («; (M + M r )e/2, where 

♦ 

M m D - TiAt - trAr. 

1 If Af -j- Af 7 can be made negative definite then 

« + M t )c/2 < -<* II r B» t (q> > 0). 

Equation (2.19) then becomes 

I r Vs-kc, ik n* +t(?T,) 

and using Schwartz's inequality we get after dividing by || « || 

jll«HS -*c»IKfl+*8T,U 
and therefore (using the fact that v(0) « u(0)) 

||?||<llLk(i. e -*wt) 

Cs 
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( 2 . 10 ) 


( 2 . 11 ) 


(2.12) 


(2.13) 


where the “constant" |j Tj Hat* maxo<T<t (j Tj(r) ||. 

If we indeed succeed in constructing M such that M ■+ M T is negative definite, with c& > 0 
independent of the sise of the matrix Af as it increases, then it follows from (2.13) that 
the norm of the error will be bounded for all t by a constant which is 0(h m ) where m is 
the spatial accuracy of the finite difference scheme (2.6). The numerical solution is then 
temporally stable. 

The rest of this section is devoted to this task for the case of m = 4, i.e, a fourth order 
accurate finite difference algorithm. 

Let the n x n differentiation matrix, D y be given by 


45 

-154 

214 

-156 

61 

-10 

10 

-15 

-4 

14 

-6 

1 

-1 

16 

-30 

16 

-1 



-1 

16 

-30 

16 

-1 



-1 

16 

-30 

16 


1 

12ft* 


-1 

16 

-30 

16 

-1 



-1 

16 

-30 

16 

-1 

1 

-6 

14 

-4 

-15 

10 

-10 

61 

—156 

214 

-154 

45 


(2.14) 

The upper two rows and the lower two rows represent non-symmetric fourth order accurate 
approximation to the second derivative without using boundary values. The internal rows 
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are symmetric and represent central differencing approximation to u** to the same order. 
Note that D is not negative definite, and neither is the symmetric part of £(£ + D T ) which 
is given by: 


90 

-144 

213 

-156 

61 

-10 

—144 

-SO 

12 

13 

-6 

1 

213 

12 

-60 

32 

-2 

0 

-15$ 

13 

32 

-60 

32 

-2 

61 

-6 

-2 

32 

-60 

32 -2 

-10 

1 

0 

-2 

32 

-60 32 


1 

24 A* 


32 -60 

32 

-2 

0 

1 

-10 

-2 32 

-60 

32 

-2 

-6 

61 

-2 

32 

-60 

32 

13 

-156 

0 

-2 

32 

-60 

12 

213 

1 

-6 

13 

12 

-30 

-144 

-10 

61 

-156 

213 

-144 

90 


(2.15) 


In order to construct M we need to specify Al, Ar , r* and tr. We construct Ai as 
follows: 


A L .Al L ' + c L AM 


(2.15) 
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where 


The o’s are given by 



O} os a* cr( 0 ... 0‘ 

°i ftj <>4 05 0 ... 0 

Q) O} Q) Q 4 Os 0 ... 0 

Oi O] Os 04 tts 0 ... 0 t 

Oi Os Os 04 Os 0 ... 0 


(2.17) 


l Oi ^2 <*3 ct 4 as 0 ... 0 J 


CL * diag {-20tt|/71, 0,...,0] 


(2.18) 



-1 5 -10 10 -5 1 0 ... 0 
— 1 5 —10 10 ••-5 1 0 ... 0 


L-i 5 -10 10 -5 1 0 ... 0 


(2.19) 


01 " 1+ i u+ s 7j+ n 7i+ s^ 

° a = " (*» + jll + 57! + gll) 

°* * + ~7l 4 2^2 + (2.20) 

04 " - (jw + \->i + pi + pi) 
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Note that gives a vector whose components ere the extrapolated value of v at x ® T* 
(i.e., VTi(t)), to fifth order accuracy; while gives a vector whose components represents 

(d*t>i /dx $ )h*. Since Cl (see 2.18) is of order unity, then 4 lv = (AW + ctAj L ))v represents 
an extrapolation of v to rr t to fifth order. 

Before using At in (2.11) or (2.6) we must define rt: 


where 


a * ^-jdiaglr,, r *' 0» • • • » o) 


( 2 . 21 ) 


n * 71/20, 

T* S (-94 - a*r,)/a, 

TS a (113-o s r,)/o, (2.22) 

r 4 = (— 56 — o«Tj)/a, 

** a (11 - o«r,)/o. 

The right boundary treatment is constructed in a similar fashion, and the formulae corre- 
sponding to (2.16) • (2.22) become: 






0 ay -4 ay_* ay.* ay-, 

0 ay -4 ay-s ay-* ay-, 

0 ay -4 ay -3 ay-* ay-, 

0 ay -4 ay -3 ay-* ay-, 

0 ay -4 ay -3 ay-* ay_, 

0 ay -4 ay -3 ay-* ay-, 

0 ay -4 ay -3 ay-* ay-, 

0 ay -4 ay -3 ay-* ay_, 


(2.23) 


(2.24) 
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C* = diag[0,0 0, -20a*/71) 

0 0 ... 0 1 -5 10 -10 5 -1 

0 0 ... 0 1 -5 10 -10 5 -1 
• 

0 0 ... 0 1 -5 10 -10 5 -1 

The a's are here: 

= 1 + 12 <y,t+ 24' ri + l2^ + 24^ 

«*-l - “( 4 ^ + T^ + l 7 * + | 7 *) 

o.v-a = 37j» + ~^ + 27 j + i^ 

/4 7 « 7 • 1 4 \ 

“• v - 3 = -(.i w+ r* + r* + 6^) 

“ j'w + 5j7j» + j'>* + jj-rA, 

^ "* •mTVM. TN'*h ^vj> 

tx = 71/2oy 



(2.25) 

(2.26) 

(2.27) 

(2.2S) 
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m-\ - (-94 - o\_it.v)/qn 

tk-i *= (113 — <xn-*tn)/<*n (2*29) 

rw-3 = (-56 - orjv-sTw)/o,y 
Tv-4 = (11 “ Q | .V-4TN , )/ Qf .V 

We «re now reedy to construct 

+ + - WAi 1 ' + e t 4‘>) + + cA***)] 

- ln.(At» + «i.A< 1 ») + »),(A«» + c*A«*>)] T } (5.30) 

Upon using equations (2.14)-(2.29) in (2.30) one gets: 


<W + A/ T 



0 

0 

-2 0 
32 -2 

-2 32 -60 32 

-2 32 -60 

-2 32 


0 


0 

-2 

32 -2 
•60 32 -2 

-2 32 -60 32 -2 0 ... 

-2 32 

-2 
0 
0 



(2.31) 
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where and are 6x6 blocks given by: 



. »J‘> + K'J 1 ' 

(2.32) 

10*) „ „**> + n^R) 

(2.33) 

i • 1 or j • 1 | 



1 < ij < 5 (2.34) 


+ QjTi) ij ? 1 J 



0 « m A* or j « JV 

”(<>>’-» T,v-j + or s-iTs-i) 


0< -V-*,A’-;<4 


( 2 . 35 ) 




-1 

0 

0 

0 

0 

O' 

0 

-30 

12 

13 

-6 

1 

0 

12 

-60 

32 

-2 

0 

0 

13 

32 

-60 

32 

-2 

0 

-6 

-2 

32 

-60 

32 

0 

1 

0 

-2 

32 

-60 J 

-60 

32 

-2 

0 

1 

O' 

32 

-60 

32 

-2 

-6 

0 

-2 

32 

-60 

32 

13 

0 

0 

-2 

32 

-60 

12 

0 

1 

-6 

13 

12 

-30 

0 

0 

0 

0 

0 

0 

-1. 


(2.36) 


(2.37) 
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The next task is to show that M « |(A#+ M*) is negative definite. We write the symmetric 
matrix A# as a sum of five symmetric matrices, 

* = + 2 * , + (24 -*>)*» + *• + **]• (2.38) 

We shall show that Afj is negative definite, and that Af,(j> = 2, ... 5) arc non-positive definite. 
The J&’ s are given by 

0 0 ] 

0—2100 
0 1-210 
0 0 1-21 

Mx * 0 0 0 1 -2 1 * + /0r, + A#f (2.39) 

• • « 

1-2 1 0 

1 -2 0 

l « ® -*J 

where Mf = | q j , Af* * | o -\/20o ] ^ Afj is the remaining (.V -2) x (A* -2) 

middle block. 
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0 


0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 -1 2 -1 

0 0 0 0 2 -5 4 -1 

-1 4-6 4 -1 

« ( 2 . 40 ) 

-1 4-6 4 -1 

-1 4 -5 2 0 0 0 0 

0 -1 2 -1 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 . 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 -1 1 

0 0 0 0 0 1 -2 1 

1 -2 1 

& * ( 2 . 41 ) 

1 -2 1 

1 -1 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 
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-26 + 0 

' 

-2 

0 

1 

28-/9 

-58 + 2/8 
- 2 or,v- 4 tfc ’_4 

32-/8 

— (®£A*-ST>f-4 
+OJV-4 ^-s) 

-2 

+<*N-4 r X-l) 

-6 

“(<*.V-iTN-4 

+«W-4Tv-») 

-2 

32-/8 

“ (<*JV-STV-4 
+OfJV-4*fc-S 

-60 + 2/9 
-2ojv-s»V-s 

32-/9 

+QW-STN-*) 

13 

-(<*X-lTs-s 
+o^_ 3 ry_ 1 ) 

0 

-2 

— (aje.ss>r-« 
+a.v-4T.v-a) 

32-/8 

-(ov_ 2 T)v_a 

+«/v-sTjv.t) 

-60 + 2/8 
-2a.v-jry.j 

12-/8 

“(flW-tTN’-S 

+a.v-sTjv.,) 

1 

-6 

-(OtN~\Tf/-A 

+<*N-4TK-\) 

13 

— (a/V-lTV-S 
+O.V-3T^-l 

12-/8 

+o.v-aTN’-i) 

-30 + 2/8 

0 

0 

1 

0 

0 

1 

t 

0 

i 


- 1/2 


Let vs consider M\ * see (2.39); it may be decomposed as follows: 


1 -1 


• • 


-1 •. 


M x = 


-1 

1 


-1 1 


The lest matrix in non-positive definite. The first term is a product of a regular matrix with 
its transpose, hence its negative is a negative definite matrix. Thus we established that M\ 



is negative definite for any finite dimension N. All its eigenvalues are negative. It remains 
to show that the eigenvalues of Mi/k* (see (2.38) are bounded away from aero by a constant 
as h -» 0 (N -+ oo). 

C o nsider a symmetric tridiagonal matrix S with, like JS fj, constant diagonals: 


b a 0 
aba 
0 aba 


(2.45) 


aba 

a b 


Designate by Dj the determinant of the upper-left j x j sub-matrix. Thus D\ » 6, Dt * 
det ‘ 

We have then D\ « b, D* * 6 a - « a and in general 

Dj * bDj-x - (2.46) 

It can be shown (see Appendix 1) that the solution to the recursion relation (2.46) is 


where 


« If A B 


* - ‘ 


(2.48) 




(2.49) 


* — [(Da — W)j)/*i -f Dj] 

Mt “Ma 


(2.30) 


B * —J— -((Da-JDjJ/ia + D,) (2.51) 

We have already shown that M\ is negative definite. The eigenvalue of jidfi are found from 

detfjfir, - IX) - (- jL _ a) . d^Ajr, _ xt) ■ (- JL - a) - 0 (2.52) 

thus either A = — 1/2& < 0 (because fo will be taken positive) or A s eigenvalue of A/i < 0. 
We would like to investigate the behavior of the eigenvalues of jjgjA/i. In particular we 
would like to show that these eigenvalues (which are negative) are bounded away from zero. 
To show this we analyse the behavior of - A/ as N increases. We now take S m Afj - A/. 

Its determinant is given by D*. a . Substituting (2.48)-(2.51 ) into (2.47) with j m N - 2 we 
get after some elementary manipulations 

= (253) 

where 


p * v4 — 6*; 6 * - 2 — A; o = 1 
r * ^1* + p* *» 2 
9 = tan-V*) 


FVora (2.52) we requite 


Dtf-a = 0 
18 


(2.54) 


(2.55) 




This is equivalent, see (2.53), to requiring 

k* 

9 ~ n - r 

Rom the definition of $ and (2.54) we obtain 


A = l,...,.V-2. 


4 / far \ V=3 ®Tm) 4 

\JV-l) 2 + A ’ ( A <0). 

Squaring (2.5?) we get a quadratic equation for A, the solution of which is 

A - -2[1±(. + W(^))' ,,J J 

For any fixed S , the smallest values of |A| is given by (2.58) for As 1. 


As A r increases, we have 


Abm * min|A| « -2 |l - cos ( y^ )] • 

“ - -2 |l - (l - 2 (A - 1)» + 0 (/?«))] 
•a 

S — II _]f^ 

(AT ~ 1)» *• 


(2.58) 


(2.59) 


(2.60) 


Thus the eigenvalues of Xfri/24 A* (and hence of Afi/24h a ) are bounded away from zero by 
the value - 

We now consider M a . One can verify that 


Af a *s 


(2.61) 


•4 




where 


0 

0 

1 0 

2 1 0 

1 -2 1 

1 -2 1 


1 0 

0 1—2 10 

0 1-2 0 
0 1 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

Therefore X/j is non-positive definite. In a similar fashion A/s is non-positive definite because 

A/s = -A/* A// (2.63) 
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with 


0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 


-1 

1 -1 


1 

1 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 


( 2 . 64 ) 


The matrices Af« and Ms are N x N 


matrices with zero entries except for 6 x 6 upper-left 


(lower-right) blocks. It is sufficient to show that these blocks are negative definite. This 


was done symbolically using the Mathematics software and plotted for 0 < 74 , 71 * < 1 and 
So * 1 . jft 4 and Ms are indeed negative definite for, 0 < 7 a, 7 t, < 1 . Thus we have shown 
that M * M + M T ) is indeed negative definite, and its eigenvalues are bounded away 
from aero by (— x a /24), even as N — * oc, and the error estimate (2.13) is valid. 


3 The Two Dimensional Case 


We consider the inhomogeneous diffusion equation, with constant coefficients, in a domain 
ft. To begin with we shall assume that ft is convex and has a boundary curve dft € C a . 
The convexity restriction is for the sake of simplicity in presenting the basic idea; it will be 
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removed later. We thus have 


lr =t (f? + 0) +/(I ’ 1 ' ;<)i *’*‘ n! ‘- 0- ' * >0 

u(x»y,0) * 

«(*>|M)I«> * *s(t) 

We shall refer to the following grid r epresentation: 


KM, 

k 

km3 

krn2 

M 


jml jmijmS j 

Figure 2: Two dimensional grid. 

We have Mn rows and M t columns inside 0. Each row and each column has a discredited 
structure as in the one 1-D case, see figure 1. Let the number of grid points in the I th row 
be denoted by A* and similarly let the number of grid points in the j* k column be C y Let 




(3.1a) 

(3.1b) 

(3.1c) 
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the solution projection be designated by £/,•,*(!). By U(t) we mean, by analogy to the 1-D 
case, 


U(<) 


(“Ui • • • » ««„i> • • • » • • • » «i M K , uajtfjt, . . . ur* h Mr) 

(ui,u 2 ,...,Uv a ) (3.2) 


Thus, we have arranged the solution projection array in vectors according to rows, starting 
from the bottom of Q. 

If we arrange this array by columns (instead of rows) we will have the following structure 
u"'(l) » (u,., , u^, . . . , ttue,; ttw, ti 2 4» • • • * «».**; • • • ? «w«,a* • • • » ) 


= 



(3.3) 


Since U^(<) is just a permutation of U(t)> there must exist an orthogonal matrix P such 
that 

U<‘>(t) . PU (3.4) 

If the length of U(t) is f, then P is an t x i matrix whose each row contains / — 1 acres and 
a single 1 somewhere. 

The second derivative operator d*(dx 7 in (3.1a) is represented on the P* row by the 
differentiation matrix i>i*\ whose structure is given by (2.14). Similarly let iP/dy* be given 
on the j* column by D\*\ whose structure is also given by (2.14). With this notation the 
Lapladan of the solution projection is: 

(•p + p) M«) » C'U + SMUM + x<*> + X</> (3.5) 
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where 


d<*' ] r ds»> i 

5*' 1 “ fl?’ ;IX»> — B}»> (3.«) 

l *>SU l Bjg J 

where DM end DM ere (t x f) matrices end have the block structures shown. TM *nd TM 
are the truncation errors associated with DM and TW, respectively. We now call attention 
to the fact that DM and DM do not operate on the same vector. This is fixed using (3.4): 

7V(f) * V*U » (D<*> + P T D r »>P)U + TM + P t T</> (3.7) 

Thus (3.1a) becomes, by analogy to (2.5), 

^ « h(D^*> + P r I>Mp)U + jk(T<*) + P T 7l») + f(t) (3.8) 

where f(t) is /(x,y;<) arranged by rows as a vector. 

Before proceeding to the semi-discrete problem let us define: 

A/j** — 0*** ~ — TRtAfit (3.9) 

where are the rj, and At defined in section 2, appropriate to the k ,tl row; similarly 

for and Ar*. In the same way, define 

Mj 9) « Dj f) - t B} A L} - rr^n, (3.10) 

where B and T stand for bottom and top. 

We can now write the semi-discrete problem by analogy to (2.6) 

~ * *(At<«> + P T M {t) P)V + *GM + kP T G M + f(t) (3.11) 
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where V is the numerical approximation to U; 

m,'*i i r *f<»> l 

X 1 ' 1 * A/i* 1 \M W m j (3.12) 

*SJ l *•&. 

and 

G (t) = + n'.prj, . • • ,(T»,ga, + Tr,gr,), .... + rrw.gTvjj • (3.13) 

Subtracting (3.11) from (3.8) we get in a fashion similar to the derivation of (2.9): 

~ = k\M M + P t M^P)E + *T, (3.14) 

where E = U — V is the two dimensional array of the errors, , arranged by rows as a 
vector. Tj is proportional to the truncation error. 

The time change of || E ||* is given by 

^ It E || a = *(E, + /> r A<‘*>/>)E) -f *(E, T 3 ) (3.15) 

The symmetric part of + P^ MW P is given by 

i|(A<<*> + + P r (M™ + M mT )P] (3.16) 

Clearly MW + AfW T and + M^ T are block*diagooal matrices with typical blocks 
given by + M^ T and AfJ rt + Af}***. We have already shown in the one dimensional 
case that each one of those blocks is negative definite and bounded away from zero by r*/24. 
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Therefore the operator (3.16) it tlto negative definite end bounded away from zero. The 
rest of the proof follows the one dis&^sdcnal case and thus the norm of the error, jj E ||, is 
bounded by a constant. 

If the domain 0 is not convex or amply connected then either rows or columns, or both, 
may be ‘interrupted* by dft. In that case the values of the solution on each “internal” 
interval (see figure (3] below) are taken as stperaie vectors. 



Decomposing “interrupted” vectors in this fashion leases the previous analysis unchanged. 
The length of U (or U<*>) is again l, where t is the number of grid nodes inside ft. The 
differentiation and permutation matrices remain ixt. Note that adding more “holes” inside 
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dU does not dung* the general approach. 


4 Numerical Example 

In this section we describe numerical results for the following problem: 

^ ~ *f Qy») 4 /(*< y, Pi (z,y)eft, t > 0, (^-1) 

where ft is the region contained between a circle of radius n> as 1/2 and inner circle of radius 
n < 0.1. The inner circle is not concentric with the outer one. Specifically ft is described by 

{(z - .5)* + (y - .5) a < 1/4} H {(* - .6)* + (y - .5)* > (.1 - *)*; 0 < S < .l} (4.2) 

The cartesian grid in which ft is embedded spans 0 < z,y < 1. We took Az = Ay, and ran 
several cases with Az *= 1/50, 1/75, 1/100. The geometry thus looks as follows: 

1 

OS 
r-0 
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The source function /(x,y, <) was chosen different from xero so that we could assign an 
exact analytic solution to (4.1). This enables one to compute the error Ea = 

“exactly" (to machine accuracy). We chose k = 1 and 

t»(x, y, t) * 1 + cos(10t - 10x 3 - 10y 3 ) (4.3) 


This leads to 


/(x, y, t) = 400(x 3 + y 3 ) cos (10* - 10x* - 10y 3 ) 

- 50sin(10t - 10x 3 - 10y 3 ) (4.4) 

From the expression for u(x,y,t) one obtains the boundary and initial conditions. 

The problem (4.1), (4.2), (4.4) was solved using both a “standard" fourth order algorithm 
(a 2-D version of (2.3)) and the new “SAT," or “bounded error,” approach described in 
Section 3. The temporal advance was via a fourth order Rungc-Kutta. 

The standard algorithm was run for Ax = 1/50 and a range of 0 <6< .01 (.09 < r< < .1). 
We found that for 6 > .0017323, the runs were stable and the error bounded for “long" times 
(10* time steps, or equivalently t *= 2). For 0 < £ < .0017233 the results began to diverge 
exponentially from the analytic solution. The "point of departure" depended on d. A 
discussion of these results is deferred to the next section. Figures 5,6,7 show the la-norm of 
the error vs. time for different radii of the inner "hole." 

The same configurations were also run using the "bounded error" algorithm described in 
Section 3 (see eq. (3.5)), and the results are shown in figures 8,9,10,11. It is seen that for 
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6'$ for which the standard methods fails, the new algorithm still has a bounded error, as 
predicted by the theory. 

To check on the order of accuracy, the “SAT* runs (with S m 0) were repeated for 
Ax * Ay m 1/75 and 1/100. Figure 12,13, and 14 show the logarithmic slope of the L}, L\ 
and !«, errors to be less than -4; i.e., we indeed have a 4* order method. That the slopes 
are larger in magnitude than 4.5 is attributed to the fact that as Ax « Ay decreases the 
percentage of “internal" points increases (the boundary points have formally only 3** oder 
accuracy). It is therefore possible that if the number of grid points was increased much 
further, the slope would tend to —4. Lack of computer resources prevented checking this 
point further. (For Ax * 0.01, running 20,000 time steps, t « .1, cpu time on a CRAY YMP 
is about 5 hours). It should also be noted that the “bounded-error” algorithm was run with 
a time step, At, twice as large as the one used in the standard scheme. At this larger At 
the standard scheme “explodes 8 immediately. 


•rr «rr 



Figure 5: 4 » 0.0017325, Standard Figure 6: 4 a 0.0017323, Standard 
scheme scheme 
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Figure 11: B * 0.0017325, SAT scheme 


Figure 12: Order of accuracy L\ 





A study of the effect of size of At shows that the instabilities exhibited above are due to 
the time-step being near the C.F.L.-lixnit. It is interesting that this C.F.L.-limit depends so 
strongly on the geom et r y . 

5 Conclusions 

(i) The theoretical results show that one has to be very careful when using an algorithm 
whose differentiation matrix, or rather its symmetric part, is not negative definite. For 
some problems, such '‘standard* schemes will give good answers (i.e., bounded errors) 
and for others instability will set in. Thus, for example, the "standard* s c h em e for 
the 1-D case has a matrix which, for all 0 < 71,71* < 1, though not negative definite 
has eigenvalues with negative real parts. This assures, in the 1-D case, the temporally 
asymptotic stability. In the 2*D case, even though each of the block sub-matrices of 
the t x i x-and-y differentiation matrices has only negative (real-part) eigenvalues, it 
is not assured that the sum of the two t x i matrices will have this property. This 
depends, among other things, on the shape of the domain and the mesh size (because 
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the mesh size determines, for a given geometry, the % end tr’s along the boundaries). 

Thus that we might have the "paradoxical* situation, that for a given domain shape, 
successive mesh refinement could lead to instability due to the occurrence of destabi- 
lizing 7*8. This cannot happen if one constructs, as was done here, a scheme whose 
differentiation matrices have symmetric parts that are negative definite. 

It is also interesting to note that if one uses explicit standard method then the allow- 
able C.F.L. may decrease extremely rapidly with change in the geometry that causes 
decrease in the 7*8. This point is brought out in figures 5 to 7 . 

(ii) Note that the construction of the 2-D algorithm, and its analysis, which were based 
on the 1-D case, can be extended in a similar (albeit more complex) fashion to higher 
dimensions. 

(iii) Also note that if the diffusion coefficient k, in the equation 

ttt = k£ 2 u 

is a function of the spatial coordinates, k * fc(x,y,x), the previous analysis goes 
through but the energy estimate for the error is now for a different, but equivalent 
norm. 
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Appendix I 

We stert with 


D, — 6Dj_i - a 7 Dj.j 


with 


Dj = 6, Dj = ft 2 - a 2 


We associ ate with (A.l) a generating function /(*), 

Multiplying (A.l) by ***** for each j > 3, and summing both sides we obtain: 

ar x 


(A.1) 


(A.2) 


(A.3) 


(A.4) 


leading to 


1 f D» i (ib ~ bD \ )x 
«* [* 2 - (t/o*)* + (!/«*). 


1 Pi + (Pa - bDi)x 

«* (* - «i t )(* - ti 8 ) 


where «i, u 2 are given by (2.48), (2.49). 


We may also present / by 


(A.5) 


/ 


1 

e a 


.(* - «,) (* - ua). 


(A.6) 
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Comparing (A.6) and (A.5) we get expression for A and B as given in (2.50), (2.51). Ex- 
panding the denominator in (A.6) we get the following series for / 



from which it immediately follows (see (A.3)) that 



(A.7) 


(A.8) 
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